home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Technotools
/
Technotools (Chestnut CD-ROM)(1993).ISO
/
lang_oth
/
extenpr
/
calpi.for
next >
Wrap
Text File
|
1988-06-16
|
3KB
|
108 lines
PROGRAM CALPI
C This program calculates pi using a method described in
C Scientific American February 1988, page 115.
C
C By Mark P. Esplin
C 16 Shawsheen Rd.
C Billerica, Ma 01821
C
INTEGER*2 Y(3000),ALP(3000)
INTEGER*2 T1(3000),T2(3000),T3(3000),T4(3000),T5(3000)
INTEGER*2 T6(3000)
INTEGER*2 NT,N
INTEGER*2 ITER,NITER,ITEMP,IFIL
CHARACTER ANS*1,FILNAM*20
DATA NT/3000/
10 WRITE(*,*) 'Enter number of desired digits'
READ(*,*) NDIG
N=NDIG/4+2
IF(N.GT.NT) THEN
WRITE(*,*) 'Not enough memory available'
GOTO 10
ENDIF
C For now number of iterations is
WRITE(*,*) 'Enter number of iterations'
READ(*,*) NITER
IFIL=0
WRITE(*,*) 'Log results to disk?'
READ(*,'(A)') ANS
IF(ANS.EQ.'y'.OR.ANS.EQ.'Y') THEN
IFIL=1
WRITE(*,*) 'File name?'
READ(*,'(A)') FILNAM
OPEN(1,FILE=FILNAM,STATUS='NEW')
ENDIF
C Set beginning values for Y and ALP.
C Y=SQRT(2)-1
ITEMP=2
CALL CVI2EX(ITEMP,T1,N)
C The square root use 4 temporaries, so it is actually using T3 - T6.
CALL EXSQRT(T1,T2,N,T3)
ITEMP=-1
CALL CVI2EX(ITEMP,T1,N)
CALL ADD(T1,T2,Y,N)
C ALP=6-4*SQRT(2)
ITEMP=4
CALL CVI2EX(ITEMP,T1,N)
CALL MUL(T1,T2,T3,N)
CALL NEG(T3,N)
ITEMP=6
CALL CVI2EX(ITEMP,T1,N)
CALL ADD(T1,T3,ALP,N)
C Main loop
DO 20 ITER=1,NITER
C T1=SQRT(SQRT(1-Y**4))
CALL MUL(Y,Y,T1,N)
CALL MUL(T1,T1,T2,N)
CALL NEG(T2,N)
ITEMP=1
CALL CVI2EX(ITEMP,T3,N)
CALL ADD(T2,T3,T1,N)
CALL EXSQRT(T1,T2,N,T3)
CALL EXSQRT(T2,T1,N,T3)
C Y=(1-T1)/(1+T1)
CALL COPY(T1,T2,N)
CALL NEG(T2,N)
ITEMP=1
CALL CVI2EX(ITEMP,T3,N)
CALL ADD(T2,T3,T4,N)
CALL ADD(T1,T3,T2,N)
CALL COPY(T4,T1,N)
C INV using T4-T6
CALL INV(T2,T3,N,T4)
CALL MUL(T1,T3,Y,N)
C T1=1+Y
ITEMP=1
CALL CVI2EX(ITEMP,T2,N)
CALL ADD(T2,Y,T1,N)
C T2=ALP*(1+Y)**4
CALL MUL(T1,T1,T2,N)
CALL MUL(T2,T2,T3,N)
CALL MUL(ALP,T3,T2,N)
C T1=-(2**(2*(ITER-1)+3))*Y*(1+Y+Y**2)
CALL MUL(Y,Y,T3,N)
CALL ADD(T1,T3,T4,N)
CALL MUL(T4,Y,T5,N)
ITEMP=2**(2*(ITER-1)+3)
CALL CVI2EX(ITEMP,T6,N)
CALL MUL(T6,T5,T1,N)
CALL NEG(T1,N)
C ALP=T2+T1
CALL ADD(T2,T1,ALP,N)
C PI=1/ALP
CALL INV(ALP,T1,N,T2)
WRITE(*,*)
WRITE(*,900) ITER
900 FORMAT(' Iteration number ',I3)
IF(IFIL.EQ.1) THEN
WRITE(1,*)
WRITE(1,900) ITER
ITEMP=1
CALL WRTEXT(ITEMP,T1,N)
ELSE
CALL PRTEXT(T1,N)
ENDIF
20 CONTINUE
STOP
END